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Summary 

The cooperative agreement covered work between August 1995 and August 1997. During 
that time. Dr. Raphael T. Haftka, Dr. John H. Garcelon and Dr. Mehmet Akgiin at the 
University of Florida worked together with Dr. Stephen J. Scotti of the Thermal Structures 
Branch of the NASA Langley Research Center. The focus of the work was efficient 
approximations of structural response and sensitivity. The effort proceeded in three 
directions as follows: 

1. Development of an approximation based on work by Kirsch, extended to efficient 
sensitivity approximations and demonstrated for structural models for the High Speed 
Civil Transport. A paper on this part was written and presented (Ref. 1). An updated 
version, Attachment A, was submitted to the AIAA Journal. 

2. Preliminary development of the adjoint method for calculating sensitivity derivatives. 
Attachment B 

3. A review of method for fast exact reanalysis, Ref. 2, was submitted to a conference, 
and a draft is included as Attachment C. 
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Approximations for structural response and design sensitivities significantly reduce the 
computational cost of structural optimization. To be useful the approximation must be 
accurate and computationally efficient. Damage tolerant design puts even more stringent 
computational demands on optimization procedures and approximations due to the large 
structural changes inherent in damage. A robust approximation is evaluated and extended to 
include design sensitivities. It is implemented in an optimization procedure and tested. It is 
shown that this approximation is accurate and efficient for damage tolerant optimum design 
of complex structural models. Lastly, a qualitative measure is developed and tested that 
estimates the computational savings using the approximation. 


Introduction 

Structural optimization often requires a large number of 
structural analyses, and for complex structures this can 
become prohibitively expensive. For this reason, it is 
co mm on to use approximations to the structural analysis 
during the optimization process to reduce the 
computational cost. There are two general categories of 
approximations according to their ranges of 
applicability in the design space. Global approximations 
are valid in large regions of the design space, while 
local approximations are only valid in the vicinity of a 
point in the design space. 1 

Damage tolerant design is defined as a design that can 
tolerate the destruction of one or more major structural 
components. In composite aircraft structures, 
conventional damage tolerant design typically results in 
a significant mass penalty. However, through optimal 
design for damage tolerance, this mass penalty can be 
significantly reduced. 2 

Damage tolerant design requires a global approximation 
because of the large changes possible in the response of 
a dama ged structure. The need for an efficient and 
accurate approximation is particularly acute when 
structures are designed for improved damage tolerance 


since a large number of damage scenarios must be 
considered to avoid designs tailored to withstand 
damage at only particular locations. In that case, each 
structural design has to be evaluated by analyzing the 
undamaged structure as well as a host of damage 
configurations under several load cases typically 
different for the undamaged and damaged cases. In 
addition, for optimization we usually need the 
derivatives of the structural response of the undamaged 
and da ma ged configurations with respect to the design 
variables. 

Local function approximations are usually based on a 
series expansion about a point in the design space. 
Storaasli and Sobieszczanski used the first-order Taylor 
series for reanalysis of structures with small structural 
changes. 3 Approximate reanalysis was done in a 
fraction of the time with displacement and stress errors 
less than 16%. Noor and Lowder found that the first- 
order Taylor series expansion technique was inaccurate 
for stiffeess changes greater than 20%. 4 They also 
considered second-order Taylor series expansions, 
which they found to be more accurate but at a 
substantial increase in computational cost. 5 Taye used a 
binomial expansion to directly approximate structural 
response. 6 The accuracy of the binomial series 
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approximation depends on the magnitude of the changes 
to the structure, making it also a local approximation. 

The reduced basis technique is a global approximation 
that expresses the structural response as a linear 
combination of known independent vectors. The 
accuracy of this technique depends on the number and 
quality of the basis vectors. The number of basis vectors 
is typically much less than the number of degrees of 
freedom. Fox and Miura used as basis vectors 
previously computed displacement vectors obtained 
from previous changes to the structure. 7 Kavanagh 
employed eigenvectors of the original structure. 8 Only 
those eigenvectors with a significant energy 
contribution were included. Noor and Lowder used 
response first derivatives for basis vectors. 4 All these 
basis vectors work well in approximating global 
changes; however, accuracy is still limited for complex 
structures, and computational efficiency may still be 
improved. 

Kirsch has developed a combined binomial-reduced 
basis approximation, which has proved to be highly 
accurate even for large structural changes in truss 
structures. 9 Similar to Noor and Lowder’ s approach, 
Kirsch uses the first several terms in the binomial series 
as basis vectors. He shows good results for changes in 
the design variables of up to 700%. 

The objectives of the present paper are: (a) to 
investigate the accuracy and efficiency of the combined 
approximation for complex wing structures under 
severe damage conditions; (b) to extend the method to 
approximations of the derivatives of damaged structural 
response with respect to design variables; (c) to include 
the approximation procedure together with a sequential 
approximate optimization method; and (d) to estimate 
the computational savings that these approximations 
provide. 

The paper begins with a section presenting the problem 
statement, followed by a summary of the combined 
approximation. Next the extension of the method to 
design derivative approximation is presented, a section 
on the optimization method follows, and two example 
problems are used to demonstrate the application to 
realistic models of wing structures. 

Problem Statement 

Given an undamaged structure’s stiffness matrix, Ko, 
and the vector of applied loads, R, the displacements, 
u 0 , can be calculated using the displacement method of 
analysis 


Ko Uq = R. 

R is assumed independent of the stiffness and damage. 
Ko will change due to damage (AK) and the damaged 
structure’s stiffness matrix is, 

K = Ko + AK. 

An exact reanalysis requires solving the modified 
system of equations for the new displacements, u new , 

K u Qew — (Ko + AK) Uqct — R. (1) 

Combined Approximations 

Local approximations are based on series expansions 
about a baseline stiffness matrix and are accurate for 
small changes in the matrix. One such local 
approximation is based on the binomial series, 

^new — Uq • Ko AK Uo + [Ko 1 AK] 2 Uq — 

[Ko' 1 AK] 3 u 0 + .... (2) 

The computational effort and the accuracy of the 
approximation depend on the number of terms retained 
in the series. 

Global approximations are valid for large changes in the 
stiffness matrix. In the reduced basis method, the 
displacement vector, u, is approximated by a linear 
combination of n linearly independent displacement 
vectors. Usually n is a small number 

u * q 0 u 0 + q, Ui + q: u : + . . . + q^ u = U B q. (3) 

Ub is & collection of the n displacement vectors. 
Substituting the approximate displacements from 
equation (3) for u new in equation (1) and pre-multiplying 
by U b t yields 


U B T KU B q = U B T R. 

(4) 

The reduced stiffness matrix is introduced as 


Kr = U b t KU b , 

(5) 

and the reduced load vector is 


Ra = U B T R, 

(6) 

So that equation (4) becomes 


Kr q - Rr. 

(7) 

The reduced stiffness matrix, Kr, is n x n in size, and 
the reduced load vector, Rr, is n x 1 in size, where n « 
m. The coefficients, q, are obtained by solving equation 
(7) and the approximate displacements, u, are then 
calculated using equation (3). 
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The success of the reduced basis method depends on 
selecting an appropriate set of linearly independent 
basis vectors, Ua. 1 The selection of basis vectors has 
been highly problem dependent; however, Kirsch uses 
the first several terms of a series as basis vectors as does 
Noor. Kirsch calls this technique Combined 
Approximations (CA). 

Kirsch has demonstrated excellent results using this 
technique on simple truss and frame problems. We have 
extended the combined approximation technique by 
including approximations for structural response 
derivatives and by applying higher order 
approximations for more complex structures. Using the 
binomial series (2), basis vectors u J5 j=l.../i-7 can be 
generated by using the following recursive relationship, 
noting that Kq is available in decomposed form and the 
first basis vector is the nominal displacements: 

u 0 — nominal displacements, 

Kq Uj = -AK Uj. lt y-1, (8) 

Using this recursive relationship, generating basis 
vectors is relatively inexpensive computationally, and 
high order approximations can easily be generated. For 
the structural models presented in this paper, generation 
of each basis vector requires 2% of the time required for 
a complete analysis. A first order approximation uses 
the first two terms of the binomial series; second and 
third order approximations use three and four terms, 
respectively. 

Approximation to Design Derivatives of 
Response With Damage 

The previous section described the approximations used 
for the response of damaged configurations. For the 
design process we need also the derivatives of the 
response with respect to changes in structural 
parameters. This section describes a process that 
approximates these derivatives. 

Consider the displacement approximation given in 
equation (3), where q are coefficients determined using 
the CA approach, and u 0 are the exact displacements of 
the nominal structure. Differentiating with respect to 
design variable x (denoted using the comma notation,) 
yields the following equation: 

u, x * U B q,* + U B ,x q = qo,x «o +. * • + qn-i,x u n -i + q 0 u*,, 

+ q^i Un.i fX (9) 

The displacement response basis vectors u 0 through 
have previously been computed. In addition,u 0tx , the 
sensitivity derivative of displacements for the 
undamaged structure, is normally computed in structural 


optimization, and can be calculated using any standard 
method such as the semi-analytical method shown in 
equation 10. 

Uo, x ~ -Ko’^AK / Ax) Uq . (10) 

The other basis vector derivatives, u 1>x through , 
are expensive to calculate and are truncated from 
equation (9) to yield the following displacement 
derivative approximation: 

u, x » b 0 u 0 +. . . + K\ u n -i + b n Uq, x = V B b, (11) 

where b are the unknown coefficients. Differentiating 
the displacement equation (1) with respect to a 

structural parameter, x, yields 

(K), x u + Ku, x = 0. (12) 

Introducing equation (3) and equation (11) into 

equation (12) gives, 

KV B b = -(K), x U B q. (13) 

Premultiplying both sides of equation (13) by V B T , 

V B T KV B b = -V B T (K), x U B q. (14) 

The reduced stiffness matrix is 

K r = V b t KV b , (15) 

and (K), x is approximated by the finite difference, 

(K), x » AK / Ax. (16) 

The coefficients, b, can be determined by solving, 

K R b = -V B T AK/AxU B q, (17) 

and the displacement derivative approximated by 

u, x = V B b. (18) 

Using u, „ the stress derivatives may be computed. It is 
important to realize that although the accuracy of the 
response approximation can be improved by using more 
basis vectors (equation 8), the accuracy of the design 
derivatives is influenced by the truncation of equation 
( 9 ). 

Optimization Procedure 

The optimization problem calls for minimiz ing the 
weight, W t subject to constraints on stresses, buckling, 
and displacements and is formulated as 

Minimize W (x), 
such that g i (x) < 0. 

where g x are the constraint functions. The optimization 
is solved by a sequential linear progr ammin g (SLP) 
procedure, developed by Scotti, which employs both 
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linear constraint approximations and derivative 
approximations . 10 In addition, the SLP procedure 
employs a penalized objective function, P> given as 

P = ^ + K£max(g,,0). 

i 

Where K is a constant set to be larger than the largest of 
all Lagrange multipliers of the constraints. This 
penalized objective ensures that when no feasible 
solution exists to the linearized problem (possibly due 
to move limits), the solution found will minimize the 
magnitude of the constraint violations. 

Move limits are incorporated to ensure the accuracy of 
the linear programming problem. To reduce the number 
of design variables, design variable linking is used, 
where the finite element model is partitioned into design 
regions and within each region, all elements are sized 
by the same design variables. Additional details of the 
optimization procedure are given in Appendix A. 



All computer simulations are done using the EAL 
analysis language. The optimization procedure 
described in Appendix A has been implemented using 
EAL. During development of the optimization 
procedure and of the approximation algorithm, EAL 
offered three key features. First, the analysis database is 
open and can be easily manipulated. Second, EAL 
provides a powerful command language with state of 
the art computational features. Third, EAL allows for 
simple inclusion of user-written FORTRAN subroutines 
that can interact with the database and be controlled 
using the command language. 

The database features allow both element stiffness and 
stress-displacement matrices to be easily manipulated. 
Elements may be grouped such that design variable 
linking is facilitated. Elements from different groups 
can be collected into sets, and the resulting set of 
elements can be manipulated as an entity . 12 For 
example, damage scenarios are described using this set 
feature. A set of elements, which may be derived from 
different groups and element types, represents that 


portion of a structure that is damaged. Manipulation of 
the set allows for those elements' stiffness contribution 
(AK) to be easily identified and modified. 

The EAL command language provides a rich set of 
linear algebra and sparse matrix operations. This 
provides a good environment for prototyping, testing, 
and implementation. D ama ge is modeled as reduced or 
eliminated structural stiffness. The stiffness of the 
structure, in sparse form, is directly manipulated to 
reflect the damage for each design condition. Likewise, 
element stress-displacement matrices are manipulated 
for damaged elements prior to constraint evaluation. 
Several FORTRAN routines are used to generate data, 
such as the local buckling constraints, that are too 
complicated, or computationally inefficient to be 
implemented using the command language. Although 
the EAL code was used for the present approximation 
development, it is important to note that implementation 
of the approximation could be done in most finite 
element codes that utilize the displacement method of 
analysis. 

EAL provides CPU utilization information for 
evaluation of computational efficiency. As 

implemented, the optimization procedure can disable 
the use of approximations for the damaged conditions. 
This allows for comparison between the approximate 
analysis and a complete reanalysis. The following 
example problems use this feature to estimate the 
accuracy and efficiency of the approximation. 



Two example problems illustrate the accuracy and 
efficiency of the method. First, a simple model of the 
High Speed Civil Transport (HSCT) shows the effects 
of the approximation order on accuracy of the 
approximations near large area wing damage. The 
second example problem is a finely meshed box beam 
model, which demonstrates how accuracy and efficiency 
of the approximations are affected as the model size 
increases. In addition, the box beam model is used in a 
full optimization using the approximation. 


4 

American Institute of Aeronautics and Astronautics 




Figure 4: Shear Stress Error vs. Damage 
Magnitude. 


Figure 4 shows the shear stress error in the skin panel 
#24. Skin panel #24 is located near the fuselage, 
adjacent to the damaged area. The improvement in 
accuracy of the third order approximation compared to 
the lower order approximations is much more dramatic 
for this element. The third order approximation with 
less than a 3% stress error is much more accurate than 
first and second order approximations. Reducing the 
stress and displacement approximation errors is 
important when considering approximating design 
derivatives of response. 


Table 1: Design Derivative Approximations 



Panel #24, shear stress | 


dx / dx 

dflog t) / dflog x) 

Central Diff. 

272749 

0.028 

Approximation 

367060 

0.037 


Rib #78. axial stress I 


do / dx 

dflog a) / dflog x) 

Central Diff. 

128192 

1.72 

Approximation 

119810 

1.61 


Panel #163, normal stress j 


da** / dx 

dflog a\*V dflog x) 

Central Diff. 

221382 

0.62 

Approximation 

242030 

0.58 


Table 1 shows a comparison of the approximate 
derivative of several typical stress components with 
respect to the top wing panel design variable. Table 1 
also shows the magnitudes of the derivatives and 
logarithmic derivatives from both the difference 
approximation and the approximation method proposed 
in this work. The logarithmic derivative represents the 
ratio of the fractional change in the function to the 
fractional change in the variable. For example, a 
logarithmic derivative (d(log a) / d(log x)) of 0.03 
indicates that a 1% change in x will lead to a 0.03% 
change in cx. When the logarithmic derivative is much 
smaller than unity, the sensitivity of the parameter is 
small and the derivative is difficult to evaluate 
accurately. As can be seen from the large differences in 
the derivative approximated by the different methods 
for panel #24 and the correspondingly low values for 
the logarithmic derivative, the accuracy of either of 


these methods may be poor. The design derivative 
approximations for rib cap #78 and panel #163 show 
much better correlation with the difference 

approximation, and both have much larger logarithmic 
derivatives than those of panel #24. Rib cap #78 and 
panel #163 are also in the vicinity of the damage. A 
third order approximation using 4 basis vectors was 
used to generate these design derivative 

approximations. 
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Figure 5: Design Regions and Top Surface Damage 
Location for the Box Beam Model 

Box Beam Test Cases 

The next model is a 180-inch by 90-inch portion of a 3- 
bay, titanium honeycomb wing box beam. This model 
consists of 5,978 nodes, 480 2-node bar elements 
representing spar caps, and 6,000 plate bending 
elements representing the titanium honeycomb. There 
are a total of 29,561 degrees of freedom in thi s model. 
Two load cases are considered that correspond to 2.5g 
and -l.Og maneuvers, respectively. There are three 
design variable regions, shown in Figure 5, dividing the 
beam evenly lengthwise for a total of 24 design 
variables. The design variables consist of the spar cap 
areas and the face sheet gauges and honeycomb core 
depths of the webs, upper surface, and lower surface. 
Each design region is optimized for minimum weight 
subject to strength and panel buckling constraints. 
These result in 12,481 constraints. All subsequent 
results are calculated at the optimum design conditions. 

First, the accuracy of the approximation was examined. 
Both the constraint and constraint sensitivities were 
computed and compared with those calculated by 
performing a complete reanalysis. A single load 
condition corresponding to a 2.5g maneuver and a 
single damage case were used. The damage considered 
was debonding of a 15-inch by 30-inch section of the 
bottom surface in the center bay 2. The approximation 
was run using 2, 3, 4, and 5 displacement response basis 
vectors generated usingequation 8. These same basis 
vectors are used for the response derivative 
approximations along with the uq,* basis vector as in 
equation 1 1. 
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Simple HSCT Test Case 

Figure 1 shows a simple model of the High Speed Civil 
Transport (HSCT). This is a relatively coarse half 
model with 193 nodes, 449 2-node bar elements, 383 
membrane elements, and 129 shear panel elements for a 
total of 533 degrees of freedom. A single load case was 
used that simulates a 2.5 g pull-up maneuver. The 
damage considered consisted of a complete loss of 
structural stiffness for the top wing skin panels from the 
connected ribs and spars at the wing’s trailing edge, 
next to the fuselage. The damage was simulated by 
removing 8 membrane elements for the top wing skin 
panels and 4 bar elements for the spar and rib caps in 
this area of the wing. A single design variable for the 
top wing panel thickness around the damage area was 
considered. The design variable is associated with 22 
triangular membrane elements that include the 8 
damaged membrane elements. 

As a global measure of the accuracy of the 
approximation, the potential energy of the approximate 
analysis is compared to that of a complete reanalysis. In 
order to track the error, the stiffness of the damaged 
elements are gradually reduced to zero by a factor, ( 1 - 
a), so that a can be viewed as a damage magnitude 
measure. First, second, and third order approximation 
errors are shown in Figure 2. Although all three 
approximations show a small global error, the accuracy 
of ^^approximations can only be truly assessed at the 
local level. 


Potential Energy Error 



Figure 2: Potential Energy Error vs. Damage 
Magnitude. 



Figure 3: Axial Stress Error vs. Damage Magnitude. 

Stress errors were also checked for elements that went 
through large changes in stress due to nearby damage. 
Figure 3 shows the axial stress error in a spar cap #3. 
Spar cap #3 is located in the vicinity of the damaged 
area (see figure 1). Both the first and second order 
approximations show reasonable accuracy; however, the 
third order approximation clearly indicates superior 
accuracy for this bar element. 
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Figure 6: Representative Stress Errors 


The accuracy of the approximation for stress is shown 
in Figure 6. The percent relative error to a complete 
reanalysis is illustrated versus the number of 
displacement response basis vectors. The average von 
Mises stress error of all panel elements is shown along 
with the stress error for two panels, which are located in 
design region two in the vicinity of the damage. These 
surface panels are considered good indicators of the 
performance of the approximations for the optimal 
design because their stress constraints are active (i.e., 
greater than -0.1). They exhibit large errors for 2 and 3 
basis vectors, but they indicate good accuracy using 4 
and 5 basis vectors. 

The relative stress error for panel 803 shows a slight 
increase when using three basis vectors compared to 
two basis vectors. This is counter-intuitive to the 
general trend of expected increasing accuracy when 
using more terms in a series approximation. The 
increasing error for this element is a local phenomenon 
and can be attributed to low quality basis vectors. When 
the fourth basis vector is included in the approximation, 
the combination of these four basis vectors adequately 
approximates the displacements in the region of panel 
803. 

The relative errors for the stress sensitivity of several 
critical web elements are shown in Figure 7. These 
elements showed the largest relative error for elements 
having logarithmic derivatives greater than 0.2 and 
active stress constraints. The sensitivity is with respect 
to the face sheet thickness of an element located in 


design region two in the vicinity of the damage. The 
average sensitivity error for elements with active stress 
constraints was 0.018, 0.013, 0.011, and 0.011 for 2, 3, 
4, and 5 displacement response basis vectors, 
respectively. Although the average errors are small and 
insensitive to the number of displacement response 
basis vectors, individual element sensitivity results are 
sensitive to the number of basis vectors (Figure 7). 



Figure 7: Representative Stress Sensitivity Errors 

The box beam is optimized subject to both the 
undamaged condition with two load cases and five 
damage conditions. The basic damage conditions 
considered upper surface debonding, lower surface 
debonding, web debonding, upper spar cap chord 
fracture, and lower spar cap chord fracture, which were 
combined with the load cases to produce ten design 
conditions. The debonding considered 15-inch by 30- 
inch sections of the web, along with top and bottom skin 
panels. Two optimization runs were made. The first 
performed exact static analyses and semi-analytical 
sensitivity calculations for each design cycle. The 
second used the CA approximations with 5 
displacement response basis vectors for all damage 
design conditions (equation 3). CA based sensitivity 
approximations are done for damage conditions using 
the 5 displacement response basis vectors along with 
u 0 , x basis vectors (equation 1 1). 

Execution times for a Silicon Graphics INDIGO-2 are 
summarized in Table 2. The load/damage condition is 
shown along with the CPU times in seconds for static 
analysis and sensitivity computations. The CPU times 
for the damage conditions are averaged for each load 


Table 2: CPU Times (seconds) for Calculation of Design Constraints and Their Sensitivities. 


Damage/Load 

Condition 

Exact Static 
Analysis 

Approx. 

Analysis 

Exact Sens. 

Approx. 

Sens. 

Ratio 

Analysis 

Ratio Sens. 

■ 

751.25 

731.71 

633.29 

629.91 

0.97 

0.99 


609.00 

98.07 

410.88 

151.65 

0.16 

0.37 

-1.0G 

Undamaged 

753.47 

727.34 

403.53 

397.85 

0.97 

0.99 

-1.0G Damage 

608.74 

98.04 

410.42 

150.61 

0.16 

0.37 
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condition. The CPU times for the CA approximation are 
compared to the CPU times for a complete reanalysis in 
this table for a single design iteration. The ratios of 
analysis times show a significant difference for the 
damage conditions, but are close to unity for each 
undamaged condition. This is due to the initial nominal 
static result required by the both the CA approximation 
and the exact analysis. The approximation shows 
significant computational savings for the damage 
conditions where reanalysis is required. It is also 
important to note that the ratio of sensitivity 
computation times do not show as great a computational 
savings as the static analysis. This is due primarily to 
the stress recovery stage of calculating the constraints 
and their sensitivities. This significant computational 
stage is identical for both the approximation and the 
exact reanalysis. Thus, the computational efficiency of 
the CA based sensitivity approximation is somewhat 
diluted. 

Figure 8 shows the CPU times for exact and 
approximate damage optimization over the course of 20 
iterations. In addition, it shows the objective function 
value for both the complete reanalysis and approximate 
reanalysis optimization runs over 20 iterations. (Note 
the expanded scale.) The box beam weight from the 
approximate optimization run was within 0.2% of the 
weight obtained using an exact reanalysis. In addition, 
no constraints were violated in either optimization run 
in the final design. 



Figure 8: CPU Times and Objective Function 
Values. 


The values of the design variables during the 
optimization were also examined. Figure 9 shows the 
lower surface spar cap areas, comparing the design 
variable values determined using a complete reanalysis 
and approximate reanalysis. Figure 10 shows the lower 
surface skin gauges. Figure 1 1 and Figure 12 show the 
upper surface skin gauges and lower spar cap areas, 
respectively. 



Upper Surface - Skin Design Variables 



Figure 11: Upper Surface Skin Gauges 
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This formula was tested by modifying the box beam 
design problem for 8 design variables, 2 load cases, and 
9 da ma ge cases, which are combined into a total of 14 
design conditions. The formula predicts a savings ratio 
of 0.341, where the actual comparison yields a ratio of 
0.324. This formula can be expected to vary for other 
types of structures. 


Conclusions 

The case studies of the simple HSCT model and the 
detailed box beam show good accuracy for the CA 
approximation. The approximation reduced the cost of 
the analysis of a damaged box model by a factor of six 
and the cost of sensitivity by a factor of three. The 
lower efficiency of sensitivity calculations is due to the 
computational cost of stress recovery, which is common 
to both the approximation and exact derivative 
calculations. Overall, CPU time for the optimization 
was reduced by about a factor of two. 

For complex structural models, using five basis vectors 
appears to provide good accuracy. Although the 
displacement response values are accurate for an 
approximation using three basis vectors, higher order 
approximations are required to yield accurate sensitivity 
values. Generation of each displacement response basis 
vector required approximately 2% of the computational 
cost of performing a single analysis, so using a high 
order approximation did not impose a significant 
penalty. 

The runtime savings estimation lets designers assess the 
computational savings that these approximations 
provide. They can get a qualitative assessment of the 
computational cost of changing the number of design 
variables and the number of damage scenarios. For most 
damage tolerant airframe designs, a large number of 
damage scenarios (related to the number of load cases 
and design variables) would be typical in order to avoid 
designs tailored only to particular damage cases. The 
CA approximation makes damage tolerant design 
feasible by allowing designers to consider more damage 
cases and load cases in optimization problems by 
reducing the computational cost of each design cycle. 
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load cases. 



Many of the design variable values show some 
oscillations near the end of the optimization run. This 
phenomenon is likely caused by linearization of the 
problem, move limit size, and errors associated with the 
sensitivity calculations. Both the approximation to 
design derivatives for the damaged structure and the 
Taylor series approximation for the undamaged 
structure (see Appendix A) contribute. The value of the 
objective function is not significantly influenced by this 
phenomenon. 

The basic premise of damage tolerant design is that a 
mass penalty is typically incurred if damage is not 
considered during optimization. When this box beam is 
optimized without considering damage, the resulting 
weight is 1 1% lighter than the damage tolerant design. 
However, when this lighter structure is damaged using 
the damage condition that generated the largest number 
of critical constraints, it can only carry 51.5% of the 
applied loading. These low values indicate the danger 
of substantial weight penalties for reinforcing the 
structure against damage after the fact. 

Using the timing results from Table 2, the following 
formula has been developed to estimate the runtime 
savings using these approximations for structures 
similar in model complexity to the box beam: 

_ lc„ (l + 0.023 dv„ )+ dc n (0. 137 + 0.009 dv n ) 

“ lc n (l + 0.023 )+ Hc n (0.8 12 + 0.023 dv „ ) 

where, R is the ratio of a single approximate design 
cycle to a single complete reanalysis design cycle; dv n is 
the number of design variables; dc n is the number of 
damage design conditions, and lc n is the total number of 
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An iteration is represented by a circuit through the 
outermost loop. Iterations are either approximate or 
exact. An approximate iteration is defined as a pass 
through the outermost loop using approximate 
sensitivity derivatives. This approximation is described 
by Scotti 10 and satisfies a second order Taylor series 
relation between the new and previous constraint 
sensitivities. An exact iteration is defmed as a pass 
through the outermost loop using exact sensitivity 
derivatives, where the semi-analytic method is used for 
displacement derivatives, which are used in a first-order 
Taylor series to determine local, perturbed constraint 
values. Constraint derivatives are then determined by 
finite differences. In the present paper, the branch to the 
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“Approximate Sensitivity Derivatives” stage was 
disabled to better illustrate the effect of the damage 
approximation. 

The optimization procedure is written primarily in the 
EAL command language with the addition of three 
external processors, which are written in FORTRAN. 11 
Processors have been developed to evaluate stress and 
local buckling constraints for a number of structural 
sections common to aircraft design. Sequential linear 
programming using the MINOS V5.4 linear 
programming routine with move limits, as described by 
Scotti 10 , performs the optimization. 

This procedure is now modified to incorporate the CA. 
instead of the exact analysis and sensitivity for damaged 
configurations. 


Figure 14: Static Analysis and Constraint 
Evaluation with Approximations 



Figure 14 shows the implementation of CA in the 
“Static Analysis and Constraint Evaluation” stage of the 
optimization procedure. Design conditions are made up 
from a load case and damage scenario. Multiple load 
cases and damage scenarios are allowed, and they may 
be combined in any way. Design conditions give a 
designer control over effectively matching critical load 
and damage situations. For design conditions with 
damage, the nominal state (u 0 ) is initially undamaged 
and is exactly solved before the effect of damage is 
approximated using CA, equations 3 - 8. Figure 15 
shows the implementation of the approximation to 
design derivatives when calculating sensitivity 
derivatives, and replaces the “Exact Sensitivity 
Derivatives” stage in figure 13. It is important to note 
that the box titled “Semi-analytic/Approx. Constraint 
Derivatives” employs the same process shown in Figure 
13 for an undamaged case. Our approximation to 
displacement derivatives, equations 15 - 18, is only 
used for damage scenarios, and constraint derivatives 
are obtained from displacement derivatives in the same 
manner as the undamaged case. 
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Introduction 

This addendum describes the work performed on the adjoint method of sensitivity 
calculation. The adjoint method is an alternative to the direct method and the finite 
difference method which are widely used. Finite difference derivatives are expensive 
computationally and prone to errors. The direct and adjoint methods are two classes of 
analytical sensitivity methods. 

The adjoint method begins with calculation of an adjoint load, which depends on the 
particular response function of interest. It requires specialized programming for each 
type of response function, and this need has delayed its widespread application. However, 
for problems with multiple load cases or multiple structural configurations, such as 
damaged versions of a structure, the adjoint method can be much more efficient than the 
direct method. This computational advantage is enhanced for derivatives of stress 
functions because these derivatives can be calculated directly rather than obtained from 
derivatives of the displacements. 

The objective of this part of the study was to lay the foundation for a demonstration of the 
advantage of the adjoint method for optimization under a variety of constraints. In 
particular, a condition was established that, when satisfied by stress and strain constraints, 
leads to a particularly easy implementation of the adjoint method. Applications with the 
von Mises stress limit, which is shown to satisfy this condition, are used as examples. 



The Direct and Adjoint Methods 


The discretized equations of equilibrium may be written in terms of the stiffness matrix 
K, the force vector f, and the displacement vector u as 


Ku=f (1) 

The recovery of the stress field a from the displacement vector may be written as 

ct =Su (2) 

When we differentiate these two equations with respect to a parameter p we get 

Kiip=f =-KpU, CTp=S p u+Su p (3) 

where subscript 'p' denotes differentiation with respect to the design variable and P is 
called the pseudo load. The direct method consists of calculating displacement and stress 
derivatives from Eq. (3). Because the matrix K is available in factored form, the solution 
for displacement derivatives is much cheaper than the solution for the displacements. 
However, obtaining stress derivatives from displacement derivatives costs about the same 
as the original stress recovery, diluting the savings associated with the direct method 
compared to finite difference calculation. The method is easily implemented in finite 
element programs using finite difference evaluation of K p . 

The adjoint method tackles directly derivatives of displacement functionals of the form 
g(u,p ). The derivative of g is given as 

g P =g,p+g, u ‘ u P ( 4 ) 
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where a comma denotes a partial derivative with respect to the following subscript and a 
dot an inner product. With a bit of algebra, this can be converted into 

where u a is the adjoint displacement field, found as a solution of the system 

Ku a = g,J (6) 

The direct method requires solutions of Eq. (3) for each design variable p, each load case, 
and each structural configuration. In contrast, the adjoint method needs to calculate 
derivatives of only potentially active constraints. The number of such constraints is of the 
order of the number of design variables. Consequently, when the number of load cases or 
structural configurations is large, the adjoint method is much more efficient than the 
direct method. Furthermore, it does not have to go through the stress recovery stage, 
translating displacement derivatives into stress derivatives. 

However, there are implementation difficulties associated with the adjoint method, which 
prevent its wide-use application. The adjoint load vector g,J requires the derivatives of 
the response function with respect to displacement components. For stress components, 
or stress functions, such as the von Mises equivalent stress, g,J depends on the stress- 
displacement relationship, and hence on the details of the finite elements used. The 
information is not always readily available, and even when it is, the adjoint load must be 
programmed differently for each different finite element, an enormous programming 
burden. Fortunately, using a continuum formulation of the adjoint method instead of a 
discretized formulation can eliminate this burden. 
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Continuum Formulation of Stress Sensitivity 


For the continuum approach, the stress function is more conveniently expressed as a 
stress functional G 


G=\g{*,p)dV 


( 7 ) 


where ct is a vector whose elements may be stresses or stress resultants. The integration 
typically represents an averaging process through a small volume of the structure. For a 
truss member a is a scalar and is the axial force in the member. For a plane-stress 
element, a is the vector of stress resultants, given by 





f. 

N = 

N, 

= t D 

£ y 




r, y 


where t is the plate thickness and D is the constitutive matrix. The sensitivity of G to 
the design variable p is obtained by differentiating Eq. (7) with respect to p to yield 


G p = \ [s, P + S.a <<r. P +°.c £ p)] dv 


( 9 ) 


where a dot indicates inner product, a subscript following a comma a partial derivative 
and that without a comma a total derivative. p=t for a plane-stress element. 
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In this formulation, it can be shown (Haftka and GQrdal, 1992, p.314) that the adjoint 
load is an initial strain given by g a so that the loading depends on the stress function but 
not on the details of the finite element. This loading can be easily applied, then, in any 
program that has initial strain loading capabilities. Under the application of the adjoint 
load, the resulting displacement field is called the adjoint displacement u a . It can be 

shown that 


G, = \{g,+«yg JdV + t' 


u 


( 10 ) 


where f is the pseudo load (Eq. 3). The integrand in Eq. (10) can be simplified under 
certain conditions. If g is a homogeneous function of degree n in its arguments, that is, if 


g(ca,cp) = c"g(a,p) 


( 11 ) 


then it can be shown that 


g. a G + Pg.p ~ n 8 


( 12 ) 


If n= 0, the constraint is unchanged (Eq. 11) when both the stresses a and the design 
variable p are scaled up or down by the same factor. Here, the case where n=0 is of 
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importance since this applies to truss members and plane-stress elements as will be 
shown. The average stress in a truss member is given by 




(13) 


where L is the length and A is the cross-sectional area. Similarly, the average von Mises 
stress in a plane-stress element is given by 



(14) 


where a is the surface area of the plate and c vM is the von Mises stress defined by 


a 


&vM — 


tG 


a 3 [^.t ~ N x N y +N 2 y + 3N 2 xy f 2 


(15) 


yp 


where a yp is the yield stress. In both Eqs. (13) and (14), the integrands are homogeneous 
to the Oth degree 

Dividing Eq. (12) by p and keeping n for the time being, 




If a is such that <j/p=c p then 


( 16 ) 
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(17) 


g.c-v.p+g.p =^S 


The condition on stress is satisfied for a truss member and a plate under in-plane loading 
when the design variables are the cross-sectional area and the plate thickness, 
respectively. For the latter element, for example, N/ =N It from Eq. (8). 

Equation (10) hence becomes 


G,-I±g 

\P J 




dv+ r-u" 


(18) 


When n= 0, the integral in Eq. (18) disappears and the sensitivity is calculated easily. 
When n is non-zero, Eq. (18) is still easier to calculate than Eq. (10); the integrand is 
simply a multiple of the function g. Thus we find that for stress constraints in truss and 
membrane elements we have an additional simplification in the calculation of the 
derivatives. 


The adjoint initial strain for a von Mises constraint is given by (Eqs. 14,15) 


/„=—!— (2N-N') (19) 

2o yp ata 

where T denotes transpose and 
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N'-[v, N, -4!fJ 


( 20 ) 


Implementation in EAL 

The adjoint approach has been implemented in EAL (Engineering Analysis Language), 
for von Mises stress constraints for plane-stress elements with thickness as the design 
variable. The initial strain (Eq. 19) is implemented via dislocations within each element 
whose stress sensitivity is sought. A runstream was written for this purpose. Simple 
plate structures having known solutions were analyzed with this runstream for 
verification. The results of this runstream were compared to those of the runstream 
written by Ed Jones of EISI and were seen to be identical. Ed Jones’ approach, however, 
was found more suitable for automating the process for an arbitrary structural mesh and 
was therefore adopted from that point on. His runstream was modified and made more 
general. 

The basic steps in sensitivity analysis with the adjoint method using EAL are the 
following: 

• Given the stress state in the structure due to the external loading, calculate the initial 
strains to be imposed and the corresponding dislocations that would induce those 
strains in regions of the structure where sensitivity is sought, namely regions where 
constraints are critical. 

• Apply the dislocations to the structure to compute the adjoint displacement vector u a . 

• Compute the pseudo load P comprising forces which, in the direct method, would be 
applied to the structure to solve for the response sensitivity u p , Eq. (3). The pseudo 
load vector consists mostly of zeros except for entries corresponding to the nodes of 
the element(s) whose thickness is being changed. 

• Use u a and P in Eq. (1 8) to compute the sensitivity. 

The pseudo load can be computed in two ways: 
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• From Eq. (3), which requires the derivative of the global stiffness matrix K. That 
derivative can be obtained from finite difference and will be filled with zeros except 
for those entries corresponding to the modified element(s). The derivative will then 
be postmultiplied with the nominal displacement u. Advantage can be taken of the 
sparse nature of K p . 

• From a continuum approach, f is the load which induces a strain field of -e/t in the 
modified element(s) with s being the strain state in the nominal structure due to the 
external loads (Haftka and Giirdal, 1992, p. 309). As such, the problem here is not 
different from imposing an initial strain that causes adjoint displacements. Hence, 
dislocations can be computed corresponding to -z!t and, in turn, proceeding one step 
further, equivalent nodal forces corresponding to these dislocations. The latter forces 
will be the pseudo loads. 

The preceding steps were implemented in a simple model of high speed civil transport 
(HSCT). The pseudo loads were computed with the continuum approach. The first and 
the second approaches were checked against each other with a simpler plate structure. 
Analysis runs were made with the HSCT model under five load cases. Currently 
investigation is carried out of how various load cases govern stress constraints relative to 
each other throughout the structure. 

A copy of the runstream that implements the adjoint method for triangular E3 1 elements 
in EAL is attached. 
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